#################################################################################### 
####################################################################################
####################################################################################
# Main Text: Figure 2
####################################################################################
####################################################################################
library(here)
library(kableExtra)
library(lmtest)
library(multiwayvcov)
library(ggplot2)
library(stringr)
####################################################################################
####################################################################################
## Read in survey data (from: 8_2020_analysis_plots.R)
dat = readRDS(here("data/dat_unrestricted.rds"))
datb = readRDS(here("data/dat_restricted.rds"))

####################################################################################
## Unrestricted

# Implement SWE of ATE: Demean all covariates, fully interact
cols = c("sex","age","education_level","Arua","Bushenyi","Ibanda","Jinja","Mbale","Mpigi","Pallisa","Mbarara","Shema",
         "village_distance_HF","village_distance_HOSP","village_distance_road","village_distance_electricity", "share2006")
dat[, paste0("c_", cols)] = lapply(dat[, cols], function(x) (x - mean(x))  )
covariates = c("treatment*c_village_distance_HF + treatment*c_village_distance_HOSP + treatment*c_village_distance_road + treatment*c_village_distance_electricity + 
                treatment*c_sex + c_age*treatment + c_education_level*treatment + treatment*c_share2006 + 
                c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment",
               "c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment")

####################################################################################
## Restricted

# Implement SWE of ATE: Demean all covariates, fully interact
cols = c("sex","age","education_level","Arua","Ibanda","Pallisa","Shema",
         "village_distance_HF","village_distance_HOSP","village_distance_road","village_distance_electricity", "share2006")
datb[, paste0("c_", cols)] = lapply(datb[, cols], function(x) (x - mean(x))  )
covariatesb = c("treatment*c_village_distance_HF + treatment*c_village_distance_HOSP + treatment*c_village_distance_road + treatment*c_village_distance_electricity + 
                treatment*c_sex + c_age*treatment + c_education_level*treatment + treatment*c_share2006 + 
                c_Arua*treatment + c_Ibanda*treatment + c_Pallisa*treatment", 
                "c_Arua*treatment + c_Ibanda*treatment + c_Pallisa*treatment")

####################################################################################
####################################################################################
## H3: Full government credit index

# Unrestricted
for (i in names(dat[,c("ca_gov_index","ca_pres_index","ca_gad_index","ca_gak_index","ca_mp_index","ca_lc_index","ca_dc_index","ca_ngo_index")])){ # Dependent Variables
  for (j in covariates) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("m",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = dat))
    # Output clustered SEs (county)
    assign(x = paste("c",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = dat),
                            cluster.vcov(lm(as.formula(model), data = dat), dat$village_id)))
  }
}

# Restricted
for (i in names(datb[,c("ca_gov_index","ca_pres_index","ca_gad_index","ca_gak_index","ca_mp_index","ca_lc_index","ca_dc_index","ca_ngo_index")])){ # Dependent Variables
  for (j in covariatesb) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("bm",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = datb))
    # Output clustered SEs (county)
    assign(x = paste("bc",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = datb),
                            cluster.vcov(lm(as.formula(model), data = datb), datb$village_id)))
  }
}

coefs = data.frame(model = grep("c\\..*t$", ls(), value = T),
                   coef = unlist(lapply(grep("c\\..*t$", ls(), value = T), function(x) get(x)[2,1])),
                   se = unlist(lapply(grep("c\\..*t$", ls(), value = T), function(x) get(x)[2,2])),
                   cov = str_extract(grep("c\\..*t$", ls(), value = T), "\\.t$|\\.c$"), 
                   n = unlist(lapply(grep("m\\..*t$", ls(), value = T), function(x) nrow(model.frame(get(x))))), stringsAsFactors=FALSE)

####################################################################################
# Create labels
coefs$dv = coefs$model
coefs$dv = gsub("^(c|bc)\\.ca_dc_index(\\.t$|\\.c$)","District\nChair", coefs$dv, perl = T)
coefs$dv = gsub("^(c|bc)\\.ca_gad_index(\\.t$|\\.c$)","District\nAgency", coefs$dv, perl = T)
coefs$dv = gsub("^(c|bc)\\.ca_gak_index(\\.t$|\\.c$)","National\nAgency", coefs$dv, perl = T)
coefs$dv = gsub("^(c|bc)\\.ca_lc_index(\\.t$|\\.c$)","Local\nCouncilor", coefs$dv, perl = T)
coefs$dv = gsub("^(c|bc)\\.ca_mp_index(\\.t$|\\.c$)","Member of\nParliament", coefs$dv, perl = T)
coefs$dv = gsub("^(c|bc)\\.ca_pres_index(\\.t$|\\.c$)","President", coefs$dv, perl = T)
coefs[coefs$n == 1477,]$cov = "Full (n=1,477)"
coefs[coefs$n == 547,]$cov = "Restricted (n=547)"

# Order factor for plotting
coefs$dv = factor(coefs$dv, levels = c("Local\nCouncilor","District\nChair","Member of\nParliament","President", "District\nAgency", "National\nAgency"))
coefs$cov = factor(coefs$cov, levels = c("Full (n=1,477)", "Restricted (n=547)"))
coefs = coefs[!is.na(coefs$dv),]

# Specify the width of your confidence intervals
interval1 = -qnorm((1-0.9)/2)  # 90% multiplier
interval2 = -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
ggplot(data = coefs, aes(x = dv)) + 
  geom_linerange(aes(x = dv, ymin = coef - se*interval1,
                     ymax = coef + se*interval1, linetype = cov), lwd = 1.5, position = position_dodge(width = 1/4)) +
  scale_linetype_manual(values=c("solid", "dotted")) +
  geom_pointrange(aes(x = dv, y = coef, ymin = coef - se*interval2,
                      ymax = coef + se*interval2, group = cov), lwd = .5, shape = 21, fill = "white", position = position_dodge(width = 1/4)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + 
  theme_bw() +
  ylab("Standard Deviations") + xlab(NULL) + #ggtitle("Effect of CHP Intervention on Political Credit") +
  theme(plot.title = element_text(size = 12, face = "bold", hjust = 0.5), text=element_text(size = 10), 
        legend.title=element_blank(), legend.position = c(0.15, 0.9),legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"))

ggsave(here("plots/figure_2.tiff"), width = 7.5, height = 5, dpi = 700)
ggsave(here("plots/credit_wpd.pdf"), width = 7.5, height = 5)


rm(list = grep("^(c|m)\\.", ls(), value = T))
rm(list = grep("^(bc|bm)\\.", ls(), value = T))